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Abstract 

Operation  of  fuel  cells  is  subject  to  inherent  uncertainty  in  various  material  and  operating  parameters,  which  causes  performance  variability  and 
impacts  the  reliability  of  the  cells.  Analysis  of  the  interactive  effects  of  parameter  uncertainty  on  the  fuel  cell  performance  is  imperative  in  a  robust 
design  endeavor.  To  this  end,  a  methodology  for  simulation  of  fuel  cell  operation  under  uncertainty  is  presented  by  considering  a  one-dimensional 
nonisothermal  description  of  the  governing  physical  phenomena.  A  sampling-based  stochastic  model  is  developed,  and  parametric  analysis  is 
presented  to  elucidate  the  effects  of  uncertainty  in  several  operating  parameters  on  the  variability  of  power  density  of  the  fuel  cell.  Robust  design 
maps  are  derived  from  the  analysis  which  provide  for  selection  of  cell  temperature  and  anode  and  cathode  pressures  as  functions  of  the  input 
parameter  uncertainty  and  target  maximum  acceptable  variability  in  the  power  density. 

©  2006  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

Power  sources  consisting  of  proton  exchange  membrane 
(PEM)  fuel  cells  offer  the  advantages  of  low  operating  temper¬ 
ature  and  pressure,  quick  start-up,  and  pollution-free  operation 
[1].  The  overall  cell  performance  is  governed  by  the  critical 
issues  of  water  and  thermal  management  [2-7].  Extensive  stud¬ 
ies  on  modeling  and  computer  simulation  of  PEM  fuel  cells 
have  been  developed  towards  improved  water  and  thermal  man¬ 
agement  through  better  understanding  the  transport  and  electro¬ 
chemical  processes.  Bernardi  and  Verbrugge  [8,9],  Springer  et 
al.  [  1 0, 1 1  ] ,  Baschuk  and  Li  [  1 2] ,  and  Rowe  and  Li  [  1 3]  developed 
one-dimensional  models  for  steady  state  operation  of  fuel  cells, 
assuming  perfect  membrane  hydration.  Amphlett  et  al.  [14,15] 
studied  the  transient  response  of  a  fuel  cell  stack  by  perform¬ 
ing  a  global  heat  and  mass  balance  analysis,  and  the  details  of 
electrochemical  phenomena  inside  the  cell  were  ignored.  Two- 
dimensional  modeling  of  transport  phenomena  in  PEM  fuel  cells 
was  presented  by  Gurau  et  al.  [16],  Um  et  al.  [17],  Wang  et  al. 
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[18],  You  and  Liu  [19]  and  You  [20],  where  two-phase  flows 
in  the  fuel  cell  systems  were  also  discussed.  Mishra  et  al.  [21] 
reported  theoretical  and  experimental  study  on  the  effects  of  dif¬ 
ferent  gas  diffusion  layer  materials  and  contact  pressure  on  the 
electrical  contact  resistance.  A  methodology  for  model-based 
design  based  on  a  one-dimensional  nonisothermal  model  was 
presented  by  Mishra  et  al.  [22],  in  which  the  optimum  operating 
and  design  parameters  were  identified  using  a  comprehensive 
parametric  analysis  on  the  various  physical  and  electrochemical 
phenomena.  Mawardi  et  al.  [23]  extended  this  analysis  to  pro¬ 
vide  an  optimization  framework  to  derive  more  general  optimum 
solutions. 

The  application  of  cell  level  models  to  predicting  the  fuel 
cell  performance  is  based  on  the  assumption  that  the  parame¬ 
ters  representing  the  physical  and  electrochemical  phenomena, 
and  the  material  properties  of  the  fuel  cell  are  deterministic. 
It  must  be  realized  that  significant  uncertainty  is  inherent  in 
such  parameters  [24,25],  arising  from  sources  such  as  oper¬ 
ating  parameter  fluctuations,  inaccuracies  in  process  control, 
empirical  determination  of  the  electrochemical  model  param¬ 
eters,  and  environmental  uncertainties.  Operating  parameters 
such  as  cell  temperature,  anode  and  cathode  pressures,  relative 
humidity,  reactants  stoichiometry,  and  dry  gas  mole  fractions 
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Nomenclature 

a 

active  area  per  unit  volume  (cm2  cm-3) 

Ci 

molar  concentration  of  species  i  (mol  cm-3) 

cp  ,i 

molar  heat  capacity  at  constant  pressure 

for 

species  i  (J  mol-1  K) 

F 

Faraday  constant  96487  (Cmol-1) 

^vap 

mass  enthalpy  of  vaporization  (J  g-1) 

k 

exchange  current  density  (A  cm-2) 

bn 

catalyst  layer  membrane  phase  current  density 

(A  cm-^) 

/ 

operating  current  density  (A  cm-2) 

k 

thermal  conductivity  (J  cm-1  K) 

kp 

hydraulic  permeability  (cm2) 

Ni 

molar  flux  of  species  i  (mol  cm-2  s) 

P 

pressure  (atm) 

Pd 

power  density  (W  cm-2) 

R 

universal  gas  constant  8.314  (J  mol-1  K) 

RH 

relative  humidity 

AS 

entropy  change  (Jg_1  K) 

T 

cell  temperature  (K) 

Wi 

molar  mass  of  species  i  (cms-1) 

X 

position  co-ordinate  along  the  cell  thickness 

Xi 

mole  fraction  of  species  i 

Greek  symbols 

Q'a 

transfer  coefficient  for  anodic  reaction 

Oi  c 

transfer  coefficient  for  cathodic  reaction 

8 

membrane  expansion  coefficient 

.mem 

ew 

volume  fraction  of  water  in  the  membrane 

? 

stoichiometry 

Od 

electro-osmotic  drag  coefficient 

Oi 

fractional  coverage  of  species  i 

K 

proton  conductivity  of  the  membrane  (W  cm- 1 

K) 

X 

membrane  hydration  coefficient  (moles 

of 

water/moles  charge  site) 

Pj 

mean  value  of  parameter  j 

P'W 

viscosity  (gem-1  s) 

1 

input  parameter  under  uncertainty 

electrical  conductivity  of  species  (Q  cm-1) 

standard  deviation  of  parameter  j 

0m 

membrane  phase  potential  (V) 

<Z>S 

solid  phase  potential  (V) 

(Pi 

mass  source  for  species  i  (gem-1) 

Subscripts 

a 

anode 

c 

cathode 

cat 

catalyst  phase 

crit 

critical 

el 

electrode  phase 

1 

liquid 

max 

maximum  value 

ref 

reference 

Superscripts 
*  optimum  value 

i  boundary  i 

eff  effective 


are  subject  to  control  fluctuations  and  thus  may  be  considered 
uncertain.  The  empirical  modeling  of  the  electrochemical  phe¬ 
nomena  may  render  parameters  such  as  the  transfer  coefficients 
in  the  Butler- Volmer  equations,  also  to  be  uncertain.  The  inter¬ 
active  effects  of  the  uncertainty  in  these  parameters  affect  both 
the  transport  and  electrochemical  phenomena  in  all  regions  of 
the  fuel  cell,  and  must  be  understood. 

A  mathematical  framework  incorporating  the  interactive 
effects  of  parameter  uncertainty  on  the  performance  of  the 
fuel  cell  is  necessary  for  a  realistic,  physics-based  simula¬ 
tion  and  robust  design.  Simulations  under  uncertainty  using 
the  framework  can,  in  turn,  be  used  to  identify  the  operating 
conditions  that  maximize  the  cell  performance  with  minimum 
variability.  To  this  end,  a  stochastic  modeling  framework  is 
developed  and  illustrated  by  considering  a  one-dimensional, 
single-phase  nonisothermal  description  of  a  PEM  fuel  cell  oper¬ 
ating  on  reformate  feed.  The  stochastic  modeling  approach 
involves  quantifying  the  uncertainty  in  the  input  parameters  in 
the  form  of  appropriate  distribution  functions,  and  propagat¬ 
ing  the  uncertainty  through  a  deterministic  model  to  construct 
the  output  variability  distributions  [26,27].  The  output  distri¬ 
butions  are,  in  turn,  used  to  obtain  reliability  or  robustness 
measures. 

In  the  present  implementation,  the  fuel  cell  material  and  oper¬ 
ating  parameters  with  uncertainty  are  represented  as  Gaussian 
probability  distributions  which  are  quantified  in  terms  of  the 
mean  and  the  variance  values.  A  sampling  method  is  used  to 
generate  combinations  of  the  input  parameters  (samples)  from 
their  respective  distributions,  and  a  one-dimensional  nonisother¬ 
mal  model  is  used  to  simulate  the  fuel  cell  operation  for  each 
sample.  The  results  of  the  simulations  are  used  to  construct 
probability  distributions  of  the  power  density.  The  output  distri¬ 
butions,  in  turn,  provide  for  extracting  the  variability  information 
needed  in  a  robust  design  endeavor.  Parametric  analysis  is  per¬ 
formed  to  elucidate  the  effects  of  uncertainty  in  the  operating  and 
design  parameters  on  the  power  density  distribution  for  several 
values  of  the  fuel  cell  temperature  and  pressures  on  the  anode 
and  cathode.  It  should  be  noted  that,  although  the  focus  of  the 
paper  is  on  the  PEM  fuel  cell,  the  stochastic  analysis  methodol¬ 
ogy  presented  herein  is  readily  applicable  to  other  types  of  fuel 
cell. 

The  paper  is  organized  as  follows:  a  deterministic  PEM 
fuel  cell  model  which  forms  the  basis  of  the  stochastic  mod¬ 
eling  is  reviewed  in  Section  2,  followed  by  a  description  of  the 
sampling-based  stochastic  modeling  framework  in  Section  3. 
Section  4  discusses  the  results  of  the  stochastic  analysis  and  the 
parametric  studies  leading  to  identification  of  robust  operating 
maps. 
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2.  Mathematical  model  of  a  PEM  fuel  cell 

A  comprehensive  physical  model  for  the  transport  and  elec¬ 
trochemical  phenomena  presented  by  Mishra  et  al.  [22]  is  used  to 
simulate  the  performance  of  a  PEM  fuel  cell  operating  on  refor¬ 
mate  feed.  The  model  considers  a  typical  PEM  fuel  cell,  where 
a  polymer  membrane  is  placed  between  an  anode  and  a  cathode 
electrode  to  form  a  membrane-electrode-assembly  (ME A).  Two 
bipolar  plates  housing  the  flow  channel  are  used  to  clamp  the 
ME  A,  as  shown  in  Fig.  E  Thin  catalyst  layers  exist  between 
each  of  electrodes  and  the  membrane,  referred  to  as  the  anode 
and  cathode  catalyst  layer,  respectively.  The  cell  is  considered 
to  be  operating  at  steady  state,  and  since  the  primary  aim  of  the 
study  is  to  present  a  framework  for  stochastic  analysis  of  fuel 
cells  using  physics-based  models,  the  discussion  is  limited  to  a 
one-dimensional  modeling  in  the  direction  along  the  cell  thick¬ 
ness.  The  modeling  further  neglects  two-phase  flow  effects,  but 
these  effects  may  be  readily  incorporated  using  approaches  such 
as  in  Refs.  [28-30].  The  modeling  includes  the  effects  of  car¬ 
bon  monoxide  (CO)  poisoning  of  the  catalysts,  as  is  prevalent  in 
fuel  cells  operating  on  a  reformate  feed.  The  governing  equations 
for  the  electrode,  catalyst  and  membrane  regions  follow  those 
given  in  Ref.  [22]  and  are  briefly  reviewed  in  the  subsections 
below. 

2.7.  Electrodes 

Fuel  cell  electrodes  are  typically  made  of  porous  carbon  paper 
or  cloth,  which  serves  to  transfer  the  reactant  species  and  to 
conduct  electrical  current.  The  mathematical  model  is  obtained 
by  considering  the  conservation  of  the  species  and  energy.  Fol- 
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Fig.  1.  Schematic  illustration  of  a  proton  exchange  membrane  (PEM)  fuel  cell. 


lowing  the  development  by  Rowe  and  Li  [13],  the  governing 
equations  may  be  written  as 


species  : 


m 

dx 


COi 


energy  : 


d  2T 

dx2 


+ 


£TOcp,  +  NiW3cpA 


d  T 
dx 


T  ^vap^3 


where  the  subscript  i  =  1,. .  .,3  denotes  the  ideal  gas  species  i,  A; 
the  molar  flux  in  the  v-direction  in  Fig.  1,  co;  the  mass  source 
term,  Wi  the  molecular  weight,  keff  the  effective  thermal  con¬ 
ductivity,  T  the  temperature,  cpj  the  specific  heat  at  constant 
pressure,  N\  the  molar  flux  of  liquid  water,  hYSL p  the  enthalpy  of 
vaporization  for  water,  I  the  current  density,  crff,  the  effective 
electrical  conductivity,  and  the  summation  JA  is  performed  with 
respect  to  all  the  gas  species  in  the  mixture. 

In  this  study,  the  feed  streams  are  considered  to  be  refor¬ 
mate  fuel  at  the  anode  side  and  humidified  air  at  the  cathode 
side.  Hence,  the  gas  species  i,  i  =  1,. .  .,3,  are  defined  as  1  =02, 
2  =  N2,  3  =  H20(g)  (i.e.,  water  vapor)  for  the  cathode  electrode 
and  1  =  H2,  2  =  CO2,  3  =  H20(g),  and  4  =  CO  for  the  anode  elec¬ 
trode.  Since  no  electrochemical  reaction  occurs  in  the  electrode 
regions,  the  mass  source  term  cc>i  is  nonzero  only  for  the  water 
vapor  species  i  =  3.  To  solve  the  water  vapor  flux  A3,  the  mass 
source  term  co 3  (corresponding  to  the  vaporization/condensation 
of  the  water  species)  is  determined  in  terms  of  the  tempera¬ 
ture  and  species  concentration  V3 ,  which,  in  turn,  is  determined 
by  the  Stefan-Maxwell  equation  [13].  Furthermore,  the  electri¬ 
cal  potential  in  the  electrode  solid,  <7>s,  is  calculated  by  Ohm’s 
law,  and  the  unknowns  in  the  electrode  regions  are  T,  A3,  X3, 
and  0S. 


2.2.  Catalyst  layers 

Catalyst  layers  are  considered  to  be  a  mixture  of  membrane, 
platinum  catalyst  (solid),  and  void  space  in  this  study.  Elec¬ 
trochemical  reactions  in  the  catalyst  regions  are  coupled  with 
the  transport  of  mass  and  energy,  resulting  in  a  potential  gra¬ 
dient  across  the  cell.  The  CO  and  H2  molecules  in  the  fuel 
compete  with  each  other  for  the  vacant  catalyst  sites,  and  high 
concentration  of  CO  may  prevent  the  electrochemical  reaction  of 
the  hydrogen  molecules,  leading  to  the  so-called  CO  poisoning 
effect. 

In  this  subsection,  the  subscript  notation  for  the  species 
involved  in  the  anode  catalyst  layer  is  1  =  H2, 2  =  H+,  3  =  H20(l), 
4  =  CO,  5  =  CO2,  and  that  for  the  cathode  catalyst  layer  is  1  =  O2, 
2  =  H+,  3  =  H20(1).  The  governing  equations  for  the  various 
electrochemical  and  transport  processes  are  derived  by  the  appli¬ 
cation  of  conservation  laws  for  the  species  and  energy,  along 
with  the  Butler- Volmer  equation  for  the  electrochemical  reac¬ 
tions,  the  Nernst-Planck  equation  for  the  flux  of  aqueous  species 
in  the  membrane,  and  the  Ohm’s  law  for  electron  transfer  in  the 
solid. 
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The  conservation  equations  of  species  and  energy  in  the  anode 
catalyst  layer  account  for  the  electrochemical  reactions  of  CO 
and  H2  as  follows: 


species  H2  : 


dAh  =  Mx) 
dx  IF 


species  H+  : 


dim 

dx 


=  j  1W  +  Mx) 


species  H20(l)  : 


dN3 

dx 


dA^4 

species  CO  :  - 

dx 


_  Mx) 

IF 

Mx) 

IF 


species  CO2  : 


dN5 

dx 


M*) 

IF 


(4) 

(5) 

(6) 
(7) 


energy  : 


d  2T 

dx2 


+ 


E 

i=  1,3-5 


NiCpjWi 


d  T 
dx 


+ 


Mx)  +  Mx) 

IF 


( TAS ) 


—  (Mx)  +  j4(x))(<Ps  —  <Pm)  +  (8) 

K 

where  the  reaction  rates,  j\(x)  and  j4(x)  for  H2  and  CO,  respec¬ 
tively,  may  be  given  by  the  Butler- Volmer  equation  as 


with  the  mean  value  of  unity,  which  reflects  both  empirical  inac¬ 
curacy  of  the  parameters  in  the  Butler- Volmer  expressions  and 
the  run- to-run  variations  in  the  operation  of  the  fuel  cell. 

Since  both  CO  and  H2  are  oxidized  to  produce  H+  [22,31], 
the  species  equation  for  proton  flux,  Eq.  (4),  includes  the  sum  of 
j i(x)  and  j4(x).  The  species  consumption/generation  in  the  elec¬ 
trochemical  reaction  of  carbon  monoxide  is  represented  by  Eqs. 
(5)-(7).  In  the  energy  equation,  Eq.  (8),  the  summation  includes 
four  non-ionic  species,  and  the  sum  j\  (x)  +  j4(x)  is  also  intro¬ 
duced  for  the  overall  anode  electrochemical  reaction  involving 
H2  and  CO.  Note  that  the  reaction  rate  j\  (x)  is  proportional  to  the 
coverage  of  hydrogen  molecules,  0\ ,  which  is  defined  as  the  frac¬ 
tion  of  the  catalyst  reactive  surface  area  covered  by  the  adsorbed 
hydrogen.  Similarly,  the  reaction  rate  of  the  CO  species,  j4(x ), 
is  considered  to  be  proportional  to  the  CO  coverage,  04,  and  CO 
concentration  c4  in  Eq.  (10).  Considering  the  adsorption,  des¬ 
orption,  and  reaction  processes  of  the  H2  and  CO  species,  the 
coverage  0\  and  04  may  be  obtained  from  a  kinetic  analysis  for 
mass  balance  [31]. 

The  six  conservation  equations  for  the  anode  catalyst  layer, 
Eqs.  (3)— (8)  introduce  four  additional  unknowns,  namely,  c  1,  c4, 
0S,  and  0m,  which,  in  turn,  are  determined  by  the  Nernst-Planck 
equation  and  Ohm’s  law  [13].  Note  that  10  unknowns  are 
involved  in  the  anode  catalyst  region:  N\ ,  im,  V3,  N4,  N5,  T, 
c  1,  c4,  <Pm,  and  @s.  The  governing  equations  for  the  cathode 
catalyst  layer  follow  those  given  by  Rowe  and  Li  [13],  and  are 
not  repeated  here  for  brevity. 


j  1  (x)  =  aifei 


ci 

^ref 

C1 


exp 


—  exp 


74  W  =  aif04 


exp 


—  exp 


(10) 


The  parameter  im  is  the  catalyst  layer  membrane  phase  cur¬ 
rent  density,  which  is  related  to  the  proton  molar  flux,  N2,  via  the 
Faraday  constant,  F,  as  im  =  FN2.  In  the  energy  equation,  Eq.  (8), 
AS  is  the  entropy  change  for  the  cathode  reaction,  @s  and  @m 
are  the  electrical  potential  in  the  catalyst  solid  phase  and  catalyst 
membrane  phase,  respectively,  and  /ceff  is  the  effective  electrical 
conductivity  of  the  membrane.  The  reaction  rates  j  \  (x)  and  j4  (x), 
depend  on  the  catalyst  reactive  surface  area  per  unit  volume,  a , 
the  reference  exchange  current  density,  at  the  reference  oxy¬ 
gen  concentration,  <4jef,  the  fuel  concentrations,  c\  and  c4,  and 
the  transfer  coefficients  a?a  and  ac.  Depending  on  the  charac¬ 
teristics  of  the  half-cell  reactions  and  the  material  properties  of 
the  catalyst  layers,  the  transfer  coefficients,  a a  and  ac ,  take  on 
distinct  values,  each  in  a  range  between  0  and  2.  In  Rowe  and 
Li  [13],  both  transfer  coefficients  are  taken  to  be  unity,  which 
yields  an  approximate  Tafel  slope  of  70  mV  decade-1.  In  this 
study,  the  parameters,  and  ac ,  are  considered  to  be  uncertain 


2.3.  Membrane 


The  membrane  of  a  PEM  fuel  cell  is  generally  made  of  a 
perfluorosulfonate  polymer,  which  acts  as  a  proton  conductor 
when  saturated  with  water.  Based  on  the  model  assumptions, 
only  two  species  are  present  in  the  membrane  region,  i.e.,  the 
liquid  water  and  the  proton.  Since  no  chemical  reaction  takes 
place  in  the  membrane,  the  species  conservation  indicates  con¬ 
stant  fluxes  for  both  species.  The  flux  of  water  is  determined  by 
the  net  effect  of  electro-osmotic  drag,  diffusion  due  to  concentra¬ 
tion  gradient,  and  convection  due  to  pressure  gradient.  The  flux 
of  protons  is  described  by  the  Nernst-Planck  equation,  which  is 
further  rearranged  to  be  in  the  form  of  the  membrane  potential 
in  this  study.  The  conservation  equations  of  energy  and  species 
in  the  membrane  region  read  [13]: 

d^T  d  i2 

energy  :  -keff—y  +  —(N\h\W\)  =  —  (11) 

dxz  dx  k 


potential  : 


d@m 

dx 


hn 

K 


+  8- 


RT 

Y 


(  3  )  —  +  -  Ni 

\l+5A/dx  k  \X J 
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H2O  flux : 


mem 

ew 


C 1 


kp 

fly 


(13) 


where  the  subscript  T  denotes  the  liquid  water  species,  h  the 
specific  enthalpy,  k  the  proton  conductivity  of  the  membrane,  8 
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the  membrane  expansion  coefficient,  X  the  membrane  hydration, 
D\  the  diffusion  coefficient  for  liquid  water  in  the  membrane, 
c™em,  the  volume  fraction  of  water  in  the  membrane,  kp  the 
hydraulic  permeability  of  the  membrane,  /rv  the  viscosity  of 
liquid  water,  and  rid  is  the  electro-osmotic  drag  coefficient.  Eqs. 
(1 1)— (13)  are  solved  for  the  three  unknowns:  T ,  <£m,  and  N\,  and 
the  readers  are  referred  to  Rowe  and  Li  [13]  for  the  values  of  the 
various  physical  properties  involved. 

The  governing  equations  formulated  in  Sections  2. 1-2.3  are 
solved  with  the  boundary  conditions  that  the  temperature,  com¬ 
position  of  the  reactant  gas  mixtures,  pressure,  and  flow  rate  in 
terms  of  stoichiometry  are  specified  at  the  inlets  of  the  anode 
and  cathode  electrodes  (i.e.,  at  the  points  a  and  f  in  Fig.  1). 
It  must  be  pointed  out  that  the  water  vapor  flux,  V3,  could  be 
calculated  by  considering  the  condensation/vaporization  pro¬ 
cesses  in  the  porous  electrode  regions.  However,  due  to  the 
difficulty  in  solving  two-phase  flow  in  the  porous  media,  the 
values  of  N 3  at  the  electrode/catalyst  interfaces,  i.e.,  points 
b  and  e,  are  set  to  be  10%  of  the  corresponding  flux  of  the 
reactant  mixture,  as  suggested  by  Rowe  and  Li  [13]  for  con¬ 
venience.  An  algorithm  developed  by  Fan  and  White  [32]  was 
implemented  to  predict  the  polarization  curve,  temperature  dis¬ 
tribution,  species  concentration  and  flux  under  various  operating 
conditions.  The  physical  properties  and  kinetic  data  adopted  in 
the  simulations  are  available  in  Refs.  [13,31].  The  polariza¬ 
tion  curve  can  be  constructed  by  solving  the  governing  sys¬ 
tem  of  equations  for  the  cell  potential,  <£>,  for  different  values 
of  the  current  density,  I.  The  power  density  of  the  fuel  cell 
can,  in  turn,  be  obtained  from  the  polarization  curve  as  the 
product  of  the  cell  voltage  and  the  corresponding  current  den¬ 
sity. 

The  deterministic  model  forms  the  basis  of  the  stochastic 
modeling  framework  discussed  in  the  next  section. 


3.  Stochastic  analysis  of  PEM  fuel  cells 

Stochastic  analysis  refers  to  numerical  simulation  based  on  a 
physical  model  where  some  of  the  parameters  are  uncertain.  The 
approach  is  that  of  quantifying  the  uncertainty  in  the  parame¬ 
ters,  propagating  the  uncertainty  through  a  deterministic  model, 
and  analyzing  the  resulting  output  parameter  distributions  for 
robustness  and  reliability  measures.  In  the  case  where  the  uncer¬ 
tainty  in  the  input  parameters  can  be  represented  as  probability 
distributions,  a  sampling  method  can  be  used  generate  stochas¬ 
tic  instances  (samples)  of  the  uncertain  parameters.  Through 
deterministic  simulations  for  each  sample,  output  parameter  dis¬ 
tributions  can  be  constructed  for  analysis. 

Fig.  2  illustrates  the  application  of  the  sampling-based 
stochastic  analysis  to  a  polymer  electrolyte  membrane  fuel  cell. 
Several  parameters  in  a  fuel  cell  exhibit  varying  levels  of  uncer¬ 
tainty  including  material  parameters  such  as  permeability,  poros¬ 
ity  of  gas  diffusion  layers,  transport  properties  (such  as  conduc¬ 
tivity,  diffusivity),  and  catalyst  loading,  among  others;  model 
uncertainty  associated  with  phenomenological  descriptions,  for 
example,  the  parameters  of  the  Butler- Volmer  equations  and 
operational  parameters.  The  present  analysis  focuses  on  mod¬ 
eling  uncertainty,  in  terms  of  the  transfer  coefficients,  aa  and 
ac ,  of  the  Butler- Volmer  equation,  as  well  as  uncertainty  in  the 
operating  parameters:  cell  temperature,  T\  anode  and  cathode 
pressures,  pa  and pc\  anode  and  cathode  relative  humidity,  RHa 
and  RHC;  anode  and  cathode  stoichiometry,  £a  and  dry  gas 
mole  fractions  at  the  cathode  and  anode,  vn2/o2  and  -Xco2/h2, 
respectively,  for  a  total  of  1 1  uncertain  parameters.  Other  types 
of  uncertainty  can  also  be  included  following  the  methodology 
described  here. 

The  uncertainty  in  each  of  the  operating  parameters  is 
described  by  a  probability  distribution  function,  and  quantified 
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Fig.  2.  Schematic  diagram  of  a  sampling-based  stochastic  analysis. 
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by  the  distribution’s  mean  value  (/r),  which  denotes  the  nomi¬ 
nal  value  of  the  uncertain  parameter,  and  its  standard  deviation 
(cr),  which  is  proportional  to  the  uncertainty  in  the  parameter. 
A  degree  of  uncertainty  can  be  expressed  in  terms  of  the  coef¬ 
ficient  of  variance  (COV),  defined  as  a!(i\  thus,  a  deterministic 
parameter  with  no  uncertainty  corresponds  to  a  COV  of  zero, 
while  a  distribution  with  a  large  uncertainty  has  a  high  coeffi¬ 
cient  of  variance.  The  choice  of  the  probability  density  function 
for  each  parameter  depends  on  the  nature  of  its  uncertainty,  and 
could  be  one  of  normal,  lognormal,  triangular,  uniform,  Weibull, 
or  other  forms  as  appropriate.  Experimental  characterization  of 
uncertainty  is  needed  to  appropriately  define  the  distribution 
function.  In  the  absence  of  such  detailed  characterization,  a 
normal  or  Gaussian  distribution,  which  is  commonly  used  to 
represent  uncertain  parameters  of  a  physical  system,  is  consid¬ 
ered  to  characterize  the  uncertain  fuel  cell  model  and  operating 
parameters,  collectively  denoted  as  §/,  i  =  1,. .  .,ra,  in  Fig.  2. 

The  combinations  of  input  parameters  are  selected  from  their 
respective  distributions  using  an  appropriate  sampling  tech¬ 
nique.  The  sampling  technique  and  the  number  of  samples  are 
critical  for  the  effectiveness  of  stochastic  model  to  accurately 
represent  the  output  distributions  corresponding  to  uncertainty 
in  the  input  parameters.  Theoretically,  an  infinite  number,  and 
practically,  a  very  large  number  of  samples  is  required  to  accu¬ 
rately  represent  an  exact  continuous  probability  distribution.  In 
turn,  the  number  of  deterministic  simulations  to  be  performed 
in  a  single  stochastic  simulation  is  equal  to  the  required  num¬ 
ber  of  samples.  Since  a  single  numerical  simulation  of  a  PEM 
fuel  cell  is  computationally  expensive,  the  computational  burden 
of  a  stochastic  simulation,  which  calls  for  multiple  determinis¬ 
tic  simulations,  is  considerably  larger.  From  a  computational 
viewpoint,  therefore,  it  is  desirable  to  minimize  the  number  of 
samples  in  a  stochastic  simulation. 

One  approach  to  alleviating  the  computational  burden  is  to 
make  the  sampling  more  effective,  thus  reducing  the  required 
number  of  samples  for  the  stochastic  analysis.  Stratified  sam¬ 
pling  methods  such  as  latin  hypercube  sampling  (LHS)  [33] 
offer  this  advantage  over  Monte  Carlo  techniques  that  are  based 
on  truly  random  sampling.  In  a  stratified  sampling  method,  if 
N  number  of  samples  is  required  from  a  one-dimensional  dis¬ 
tribution,  the  distribution  is  divided  into  N  intervals  (strata)  of 
equal  probability,  and  one  sample  is  picked  randomly  from  each 
interval  to  generate  the  samples.  LHS -generated  samples  better 
represent  the  entire  distribution  compared  to  the  Monte  Carlo 
technique  in  which  the  samples  are  selected  randomly  and  may 
not  cover  the  entire  distribution.  Note  that  other  stratified  sam¬ 
pling  methods  such  as  the  Hammersley  sequence  sampling  [34] 
may  also  be  used  following  the  methodology  presented  in  this 
paper. 

The  LHS  method  is  used  to  generate  a  set  of  N  samples  from 
the  distributions  of  the  m  uncertain  parameters,  as  depicted  in 
Fig.  2.  The  selection  of  the  number  of  samples  for  this  study  is 
based  on  a  stochastic  convergence  analysis  discussed  in  the  next 
section  on  presentation  of  results.  As  noted  earlier,  the  stochastic 
analysis  presented  here  considers  1 1  uncertain  parameters,  i.e., 
m=  11.  Each  sample,  j,  represents  a  combination  of  the  uncertain 
input  parameter  values  {^jp  /=!,..  .,ra},  and  the  deterministic 


PEM  fuel  cell  model  is  used  to  simulate  the  performance  for  each 
such  sample,  j=  1,. .  .,N.  The  cell  polarization  curve  expressed 
as  the  power  density  variation  with  the  current  density,  Pd  CO,  as 
obtained  from  the  simulations  is  used  as  the  output  parameter 
for  the  stochastic  analysis. 

Based  on  the  multiple  simulations  corresponding  to  the  N 
samples,  a  distribution  of  the  power  density  variation  is  obtained 
as  illustrated  by  the  shaded  band  in  Fig.  2.  At  any  given  current 
density,  a  probability  distribution  function  can  be  constructed 
for  the  power  density,  which  is  used  to  extract  the  mean  (/xpd), 
the  standard  deviation  (a/>d),  and  the  coefficient  of  variance 
[(cr//x)p  ].  The  results  of  the  stochastic  analysis  are  presented 
in  terms  of  these  parameters  in  the  next  section. 

4.  Results  and  discussion 

Since  the  PEM  fuel  cell  model  described  in  Section  2  forms 
the  basis  of  the  stochastic  analysis,  the  model  is  validated  first 
through  comparison  with  numerical  and  experimental  data  in  the 
literature.  The  results  of  the  validation  on  the  polarization  curve 
are  presented  in  Fig.  3,  in  which  the  solid  lines  correspond  to  the 
present  simulation  results,  the  discrete  diamond  markers  denote 
numerical  prediction  by  Rowe  and  Li  [13],  and  the  circles  rep¬ 
resent  the  experimental  data  from  Springer  et  al.  [10].  The  cell 
potential  decreases  monotonically  with  increasing  current  den¬ 
sity,  as  expected,  due  to  increased  ohmic  loss  in  the  membrane. 
The  present  model  predictions  show  good  agreement  with  the 
numerical  prediction  by  Rowe  and  Li  [13]  over  the  entire  range 
of  current  density  in  the  plot. 

The  model  is  also  validated  with  the  experimental  data  from 
Springer  et  al.  [10]  where  carbon  monoxide  (CO)  poisoning  at 
100  ppm  in  the  feed  hydrogen  stream  at  the  anode  was  con¬ 
sidered.  The  cell  potential  decreases  sharply  in  the  presence 
of  CO  poisoning  in  the  anode  feed  stream,  due  to  the  adsorp¬ 
tion  of  CO  species  on  the  catalyst  surface,  which  impedes  the 
electrochemical  reaction  of  the  hydrogen  fuel.  Consequently,  a 
smaller  cell  potential  is  seen  for  a  fixed  value  of  current  den¬ 
sity.  It  should  be  noted  that  a  sharp  drop  of  cell  potential  0  is 


Fig.  3 .  Validation  of  the  numerical  PEM  fuel  cell  model  with  the  data  from  Refs. 
[13,10]. 
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Number  of  Samples,  N 

Fig.  4.  Stochastic  convergence  analysis  to  determine  the  minimum  number  of 
samples:  (a)  convergence  of  mean  power  density  with  the  number  of  samples 
and  (b)  convergence  of  standard  deviation  of  the  power  density  with  the  number 
of  samples. 

observed  between  1=  0.2  and  0.4  A  cm-2  for  the  case  with  CO 
poisoning.  As  current  density  increases,  the  loss  due  to  the  limi¬ 
tation  of  mass  transport  rate  becomes  dominant,  and  with  further 
increase  in  the  current  density,  the  partial  pressure  of  oxygen  at 
the  catalyst  layer/cathode  interface  rapidly  decreases,  preventing 
oxygen  from  reaching  the  reaction  site  and  thus  rapidly  reduc¬ 
ing  the  voltage.  The  present  model  predictions  are  seen  to  be  in 
close  agreement  with  the  experimental  data  for  all  current  den¬ 
sity  values.  The  comparisons  in  Fig.  3,  therefore,  demonstrate 
the  accuracy  of  the  deterministic  process  model,  which  serves 
as  the  basis  of  the  sampling-based  stochastic  analysis  results 
discussed  in  the  remainder  of  this  section. 

The  fidelity  of  the  stochastic  analysis  is  directly  linked  to 
the  accuracy  of  the  shape  and  moments  of  the  power  density 
distribution  obtained  from  the  simulations.  While,  theoretically, 
an  infinite  number  of  samples  is  needed  for  an  accurate  calcu¬ 
lation  of  the  moments,  faster  computation  calls  for  minimizing 
the  number  of  samples.  The  minimum  number  of  samples  for 
the  stochastic  simulations  was  determined  using  a  stochastic 
convergence  analysis  of  the  power  density  distribution.  Fig.  4 
presents  the  variations  of  the  mean  (Fig.  4(a))  and  the  standard 
deviation  (Fig.  4(b))  of  the  power  density  as  a  function  of  the 


Table  1 


Values  of  the  operating  parameters  used  in  the  studies 


Operating  parameters 

Figs.  4-6 

Figs.  7-12 

Current  density,  I  (A  cm-2) 

0.7 

Variable 

Cell  temperature,  T  (°C) 

100 

100 

Anode  pressure,  pa  (atm) 

6.5 

10 

Cathode  pressure,  pc  (atm) 

6.5 

15 

Anode  stoichiometry,  £a 

4.0 

1.1 

Cathode  stoichiometry,  £c 

4.0 

7.0 

Anode  relative  humidity,  RHa 

0.6 

1.1 

Cathode  relative  humidity,  RHC 

0.6 

0.1 

Anode  dry  gas  mole  fraction,  vco2/H9 

0.0 

0.0 

Cathode  dry  gas  mole  fraction,  vn2/Oo 

3.76 

0.0 

number  of  samples  used  in  the  stochastic  simulation.  The  results 
correspond  to  a  specific  current  density  of  0.7  A  cm-2,  and  the 
other  operating  parameters  listed  in  Table  1 .  It  is  observed  that 
both  the  mean  and  the  standard  deviation  values  converge  as  the 
number  of  samples  increases,  and,  further,  that  at  100  samples, 
the  mean  value  converges  to  within  0.1%  (Fig.  4(a))  and  the 
standard  deviation  value  converges  to  within  3.0%  (Fig.  4(b)) 
of  their  respective  values  at  larger  number  of  samples.  It  was 
found  that  the  convergence  characteristic  -  i.e.,  the  number  of 
samples  beyond  which  the  standard  deviation  (mean)  converges 
to  within  3.0%  (0.1%)  -  remained  practically  invariant  to  the 
change  in  the  initial  random  seed  used  in  the  generation  of  the 
samples  using  LHS.  Based  on  the  stochastic  convergence  analy¬ 
sis,  a  sample  size  of  100  was  selected  for  the  simulations  reported 
in  this  section. 

Fig.  5  shows  a  distribution  of  the  power  density  from  a 
stochastic  simulation  corresponding  to  a  current  density  of 
0.7  A  cm-2.  Recall  that  the  uncertain  parameters  in  the  stochas¬ 
tic  simulations  are  represented  as  symmetric  Gaussian  distri¬ 
butions;  however,  the  histogram  of  the  power  density  reveals 
a  distinct  skew  in  the  distribution,  which  reflects  the  nonlin¬ 
earity  of  the  relationship  between  the  power  density  and  the 
uncertain  parameters.  Furthermore,  it  is  important  to  note  that 
as  a  result  of  the  nonlinearities,  the  extreme  values  of  the  power 


Fig.  5.  Power  density  distribution  obtained  from  stochastic  analysis  using  100 
samples. 
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Fig.  6.  Coefficient  of  variance  of  power  density  as  a  function  of  input  uncertainty  for  various  values  operating  parameters:  (a)  cell  temperature,  (b)  anode  and  cathode 
pressures,  (c)  anode  and  cathode  stoichiometry,  (d)  relative  humidity  in  the  anode  and  cathode,  (e)  N2/O2  mole  fraction,  and  (f)  CO2/H2  mole  fraction. 


density  do  not  necessarily  correspond  to  the  extreme  values  of 
the  uncertain  input  parameter  distributions;  thus,  a  sampling- 
based  stochastic  analysis  is  better  suited  for  robust  design  than 
other  techniques  such  as  interval  analysis  that  consider  only  the 
extreme  values  of  the  uncertain  parameter  ranges.  The  bold- 
dashed-line  in  Fig.  5  represents  the  mean  value  of  the  dis¬ 
tribution  (fipd  =  0.345  W  cm-2)  and  the  dotted  lines  at  0.308 
and  0.382  W  cm-2  indicate  fipd  —  crpd  and  fipd+  crpd  values, 
respectively.  It  is  seen  that  about  75%  of  the  distribution  falls 
within  this  one- standard  deviation  band  around  the  mean.  The 
figure  suggests  that  the  mean,  /zpd,  the  standard  deviation,  crpd, 
and  the  coefficient  of  variance,  (o//x) P  ,  are  appropriate  mea¬ 
sures  of  the  power  density  distribution  and  will  be  used  in  the 
analysis  presented  here. 

Fig.  6  shows  the  coefficient  of  variance  (COV)  of  the  power 
density  as  a  function  of  input  uncertainty,  expressed  as  the  COV 
of  the  input  parameters,  (<r//x)^,  for  various  values  of  the  operat¬ 
ing  parameters  under  uncertainty.  It  should  be  noted  that  as  each 
operating  parameter  shown  in  Fig.  6  is  varied,  the  other  param¬ 
eters  assume  the  values  listed  in  Table  1.  The  results  share  a 
common  trend  that  as  the  COV  of  the  input  parameters  increases, 
the  COV  of  power  density  also  increases.  Fig.  6(a)  presents  the 
effects  of  the  cell  temperature,  T,  in  which  it  is  seen  that  higher 


cell  temperature  results  in  lower  COV  of  the  power  density  for 
all  values  of  input  COV,  (a//x)^.  It  was  reported  in  Ref.  [22] 
that  power  density  increases  with  increasing  cell  temperature. 
Thus  the  results  indicate  that  as  the  cell  temperature  increases, 
the  increase  in  the  standard  deviation  of  power  density  is  less 
than  the  increase  in  the  mean  power  density,  causing  the  COV  of 
power  density,  (cr//z)p  ,  to  decrease  as  T  increases.  The  effects 
of  varying  electrode  pressures,  pa  and  pc ,  on  the  COV  of  power 
density  are  presented  in  Fig.  6(b),  which  shows  that  lower  elec¬ 
trode  pressures  cause  smaller  power  density  variance.  However, 
it  should  be  noted  that  small  values  of  electrode  pressures  also 
result  in  small  power  density  [22],  which  may  not  be  a  desirable 
performance.  Fig.  6(c)-(f)  show  that  the  effects  of  the  varia¬ 
tion  in  the  anode  and  cathode  stoichiometry,  £a  and  relative 
humidity  of  the  anode  and  cathode,  RHa  and  RHC,  and  dry  gas 
mole  fractions,  vn2/o2  and  -xco2/h2>  on  the  power  density  vari¬ 
ance  are  less  significant  than  the  effects  of  T,  pa  and  pc,  for  all 
values  of  (cr//x)|.  For  the  remainder  of  the  parametric  study  of  the 
stochastic  analysis,  therefore,  the  paper  focuses  on  the  effects  of 
operating  cell  temperature,  T,  anode  pressure,  pa,  and  cathode 
pressure,  pc,  on  the  variation  of  power  density. 

Fig.  7(a)-(i)  present  the  mean  power  density,  /zpd,  as  a  func¬ 
tion  of  the  current  density  for  a  range  of  input  parameter  uncer- 
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Fig.  7.  Mean  and  standard  deviation  values  of  power  density  as  a  function  of  current  density  at  various  values  of  cell  temperature  and  level  of  input  uncertainty. 
Other  operating  parameters  are  given  in  Table  1 . 


tainty,  (cr//x)^  and  cell  temperature,  T.  The  variability  in  the 
power  density  as  a  function  of  the  current  density  is  illustrated 
by  the  dashed  lines,  which  represent  the  values  of  fipd  —  opd 
and  fjipd  +  apd.  For  a  given  current  density,  the  standard  devia¬ 
tion  of  the  power  density,  opd ,  is  the  magnitude  of  the  difference 
between  either  of  the  two  dashed  lines  and  the  solid  line  in  the 
plots.  The  results  show  that  as  the  current  density  increases,  the 
variability  of  power  density  also  increases  as  seen  in  the  widen¬ 
ing  of  the  gap  between  the  dashed  lines  in  the  plot.  The  increase 
is  more  pronounced  for  cases  with  higher  input  uncertainty,  as 
observed  in  Fig.  7(c),  (f),  and  (i)  (the  right  column  of  plots).  As 
(< al pc )%  increases  (from  left  to  right  in  Fig.  7),  the  peak  mean 
power  density  slightly  decreases,  while  the  standard  deviation 
of  the  power  density,  crpd ,  clearly  increases,  which  signifies  that 
increase  in  the  input  parameter  uncertainty  leads  to  a  significant 
performance  variability  of  the  fuel  cell.  As  the  cell  temperature 
increases  (from  top  to  bottom  in  Fig.  7),  both  the  mean  and  the 
standard  deviation  of  the  power  density  increase.  However,  the 
increase  in  the  mean  value  dominates  the  increase  in  the  standard 
deviation,  resulting  in  a  decrease  in  the  coefficient  of  variance, 
(<T/fjL)p ,  ,  with  temperature,  T,  as  noted  in  Fig.  6(a).  These  results 
suggest  that  the  fuel  cell  should  be  operated  at  high  temperature 
from  the  viewpoint  of  increasing  the  mean  power  density,  fi  pd ; 


however,  minimizing  variance,  opd ,  calls  for  a  low  temperature 
operation.  Identification  of  robust  operating  regimes  is  based  on 
balancing  these  two  competing  considerations. 

It  is  instructive  to  examine  the  variation  of  the  mean  values  of 
the  power  density  with  respect  to  their  corresponding  standard 
deviations  so  that  the  maximum  achievable  mean  power  den¬ 
sity  can  be  identified  for  an  allowable  maximum  variance  in  the 
power  density.  To  this  end,  Fig.  8(a)-(c)  present  the  data  from 
Fig.  7  in  this  format,  for  crpd  in  the  range  of  0-0. 10  W  cm-2, 
and  for  the  three  input  uncertainty  levels  and  the  three  temper¬ 
atures  considered  in  Fig.  7.  For  a  given  value  of  opd ,  the  mean 
power  density,  fipd,  increases  with  cell  temperature,  owing  to 
the  reduced  activation  and  concentration  losses,  and  decreases 
as  (<t//x)^,  increases  (Fig.  8(a)-(c)).  It  is  seen  in  Fig.  8(a)-(c)  that 
for  each  combination  of  temperature  and  input  uncertainty,  the 
mean  power  density  is  maximized,  /jl*p  ,  for  an  optimum  value 
of  the  standard  deviation,  Op  »  as  identified  for  an  example  case 
in  Fig.  8(a).  This  suggests  that  in  order  to  operate  the  fuel  cell  at 
its  maximum  power  density,  a  certain  degree  of  variability  rep¬ 
resented  by  <Tp  is  inevitable,  and  that  the  variability  can  not  be 
reduced  (or  increased)  further  without  also  compromising  on  the 
mean  power  density.  It  is  interesting  to  note  that  the  mean  power 
density  is  reduced  for  both  a  tighter  and  a  relaxed  specification  on 
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Fig.  8.  Mean  value  of  power  density  as  a  function  of  the  standard  deviation  for 
various  values  of  cell  temperature  and  level  of  input  uncertainty.  Other  operating 
parameters  are  given  in  Table  1 . 

the  allowable  variance  relative  to  crp ,  .  Fig.  8  further  shows  that 
the  optimum  variance  increases  with  temperature  at  all  levels  of 
input  uncertainty.  For  a  given  temperature,  crp  ,  increases  with 
(< a l/i )£,  which  indicates  that  as  the  input  uncertainty  increases, 
the  fuel  cell  can  not  be  operated  at  the  maximum  power  density 
without  corresponding  increase  in  the  variability. 

Following  the  format  in  Fig.  8,  Fig.  9  illustrates  the  effects  of 
anode  pressure,  pa  on  the  mean  and  standard  deviation  of  power 
density,  at  a  cell  temperature  of  T=  100  °C  and  cathode  pressure 
of  pc  =  15  atm.  It  is  also  observed  that  for  each  anode  pressure, 
there  exists  an  optimum  standard  deviation  of  power  density, 
crp  ,  which  corresponds  to  a  peak  value  of  mean  power  den¬ 
sity,  /Zp  .  As  input  COY,  (cr//z)£,  increases  (Fig.  9(a)-(c)),  <7p  , 


Fig.  9.  Mean  value  of  power  density  as  a  function  of  the  standard  deviation  for 
various  values  of  anode  pressure  and  level  of  input  uncertainty.  Other  operating 
parameters  are  given  in  Table  1 . 

increases,  and  the  corresponding  peak  mean  power  density,  /Zpd , 
decreases,  for  all  values  of  anode  pressure  considered.  Both  peak 
mean  power  density,  /Zp  ,  and  the  optimum  standard  deviation, 
<jpd,  exhibit  a  non-monotonic  trend  with  respect  to  the  anode 
pressure,  with  pa  =  6  atm  being  an  optimum  within  the  range 
of  anode  pressures  considered.  At  the  low  anode  pressures,  an 
increase  in pa  enhances  the  power  density  due  to  reduced  concen¬ 
tration  loss.  However,  further  increase  in  pa  causes  membrane 
dehydration  and  lower  performance  due  to  a  decrease  in  the 
water  diffusivity  in  the  anode  region  [23] .  It  should  also  be  noted 
that  the  effects  of  anode  pressure  on  the  variation  of  power  den¬ 
sity  are  not  as  significant  as  those  of  the  cell  temperature  seen 
in  Fig.  8. 
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Fig.  10.  Mean  value  of  power  density  as  a  function  of  the  standard  deviation  for 
various  values  of  cathode  pressure  and  level  of  input  uncertainty.  Other  operating 
parameters  are  given  in  Table  1 . 

The  effects  of  cathode  pressure  on  the  mean  power  density 
variation  with  the  standard  deviation  are  shown  in  Fig.  10,  for 
input  COV  of  0.02  (Fig.  10(a)),  0.05  (Fig.  10(b)),  and  0.10 
(Fig.  10(c)),  with  the  cell  temperature  fixed  at  T-  100  °C  and 
anode  pressure  at  pc  =  10  atm.  It  is  observed  that  as  the  cathode 
pressure,  pc ,  increases,  the  mean  power  density,  fipd,  increases 
monotonically,  for  any  value  of  the  standard  deviation,  due  to  the 
reduced  concentration  loss  at  the  high  cathode  pressure.  Further, 
as  previously  observed  in  Figs.  8  and  9,  the  optimum  standard 
deviation  of  power  density,  cr|>d,  also  increases  as  the  (<r//x)^ 
increases.  For  the  values  of  cell  temperature,  cathode  pressure, 
and  input  COV  considered,  both  the  peak  mean  power  density, 
/Zp  ,  and  the  optimum  standard  deviation,  Op  ,  increase  mono¬ 


tonically  with  pc.  The  results  in  Fig.  10  also  indicate  that  the 
effects  of  the  cathode  pressure,  pc ,  are  more  significant  than 
those  of  the  anode  pressure. 

From  the  foregoing  results,  it  is  evident  that  conditions  that 
lead  to  maximizing  the  mean  power  density  also  cause  the 
variability  to  increase.  Thus,  it  is  of  interest  to  determine  the  opti¬ 
mum  cell  temperature,  and  anode  and  cathode  pressures,  which 
maximize  the  mean  power  density  subject  to  constraint  on  lim¬ 
iting  the  standard  deviation  of  power  density  to  within  a  desired 
value,  for  given  input  uncertainty.  The  results  in  Figs.  8-10  can 
be  used  to  construct  design  charts  that  provide  for  determining 
the  maximum  achievable  mean  power  density  as  a  function  of 
input  COV  and  desired  maximum  standard  deviation  of  power 
density.  The  design  charts  and  the  illustration  of  the  design 
methodology  are  the  focus  of  the  remainder  of  the  discussion. 

Fig.  11  presents  contour  maps  of  the  maximum  realizable 
mean  power  density,  p />d?max,  as  a  function  of  the  standard  devi¬ 
ation,  cj pd ,  and  input  coefficient  of  variance,  (<r//z)^,  for  the  three 
values  of  the  cell  temperature.  For  all  the  three  cell  temperatures, 
the  contours  exhibit  similar  trends  with  different  magnitude  of 
power  density  values.  The  maximum  realizable  mean  power 
density  varies  from  0.56  W  cm-2  for  T=6 0°C  (Fig.  11(a))  to 
0.93  W  cm-2  (Fig.  11(c)).  The  maximum  power  density  val¬ 
ues  are  seen  to  concentrate  along  a  band  extending  from  the 
lower  left  corner  to  the  top  mid  section  of  each  plot  frame.  The 
design  plots  reveal  interesting  trends:  for  stringent  requirements 
on  the  maximum  variability,  for  example  crpd  =  0.02  W cm-2, 
the  maximum  achievable  mean  power  density  decreases  mono¬ 
tonically  as  ( alp )|  increases.  However,  if  a  greater  variability 
were  tolerable,  for  example,  (a//z)^  =  0.05  W  cm-2,  /xpd,max, 
increases  for  0  <  (cr//z)|  <  0.05  and  decreases  for  (cr//z)|  >  0.05. 

The  contours  of  maximum  mean  power  density  correspond¬ 
ing  to  different  anode  and  cathode  pressures  are  presented  in 
Fig.  12(a)-(c)  and  (d)-(f),  respectively,  for  a  fixed  cell  tem¬ 
perature  of  100  °C.  The  shapes  of  the  contours  in  Fig.  12  are 
similar  in  shape  to  those  in  Fig.  11,  in  that  the  larger  values  of 
the  mean  power  density  are  clustered  around  regions  in  each 
plot.  The  magnitude  of  the  maximum  achievable  mean  power 
density,  /zpd?max,  increases  as  the  anode  pressure  increases  from 
pa  =  2  atm  (Fig.  12(a))  to  pa  =  6atm  (Fig.  12(b)),  and  slightly 
decreases  with  further  increase  of  pa  to  10  atm  (Fig.  12(c)). 
On  the  other  hand,  the  results  for  different  values  of  cathode 
pressure  in  Fig.  12(d)-(f)  show  that  the  maximum  mean  power 
density  increases  monotonically  with  cathode  pressure,  reach¬ 
ing  a  value  of  0.93  W  cm-2  for  pc  =  15  atm  (Fig.  12(f)).  Once 
again,  the  maximum  achievable  mean  power  density  is  seen  to 
decrease  monotonically  with  (<r//z)^  at  small  values  of  allow¬ 
able,  c>pd ,  whereas  at  larger  allowable  values,  a  non-monotonic 
trend  is  evident  in  the  variation  of  p pd,max  with  (<r//z)^. 

The  use  of  the  contour  maps  in  Figs.  1 1  and  12  for  a  robust 
design  of  the  fuel  cell  operating  parameters  is  illustrated  by  con¬ 
sidering  three  design  cases — Case  A  in  which  the  variance  in 
power  density,  crpd ,  is  to  be  limited  to  0.02  W  cm-2  for  a  cell 
with  an  operational  parameter  uncertainty,  (<r//z)^  =  0.10;  Case 
B  in  which  the  variance  in  power  density,  c>pd,  is  to  be  limited 
to  0.02  W  cm-2  for  a  cell  with  an  operational  parameter  uncer¬ 
tainty,  ( alp )|  =  0.05  and  Case  C  in  which  the  variance  in  power 
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Fig.  1 1 .  Contour  maps  representing  the  maximum  mean  value  of  power  density, 
in  a  two-dimensional  field  of  the  standard  deviation  of  power  density  and  input 
uncertainty,  for  various  values  of  cell  temperature.  Other  operating  parameters 
are  given  in  Table  1 . 

density,  crpd ,  is  to  be  limited  to  0.05  W  cm-2  for  a  cell  with  an 
operational  parameter  uncertainty,  =  0.05. 

•  Corresponding  to  Case  A  specifications,  it  is  seen  from  Fig.  1 1 
that  the  maximum  realizable  mean  power  density,  /x/>d?max, 
is  highest  at  0.40  W  cm-2  for  an  operating  temperature  of 
100  °C,  and  anode  and  cathode  pressures  of  10  and  15  atm, 
respectively  (Fig.  11(c)).  Similarly,  Fig.  12(a)-(c)  reveal  that 
the  design  specifications  for  Case  A  correspond  to  a  maxi¬ 
mum  achievable  mean  power  density  of  0.40  W  cm-2  for  all 
anode  pressures  in  the  range  2-10  atm,  a  cell  temperature  of 
100  °C  and  a  cathode  pressure  of  15  atm,  and  Fig.  12(d)— (f) 
indicate  that  the  maximum  realizable  mean  power  density 


is  0.40  W  cm-2  for  pc  =  15  atm,  T-  100  °C  and  pd  =  1 0  atm. 
Based  on  these,  the  robust  operating  conditions  for  Case  A 
are  deduced  as  T-  100  °C,  pd  =  2  atm,  and  pc  =  15  atm.  Note 
that  the  solution  is  invariant  to  the  anode  pressure,  and  the  low¬ 
est  value  of  2  atm  is  chosen  based  on  minimizing  the  pumping 
requirements.  The  robust  solution  yields  a  mean  power  den¬ 
sity  of  0.40  W  cm-2  with  a  variability  of  0.02  W  cm-2,  which 
represents  a  coefficient  of  variance  of  5%  on  the  power  den¬ 
sity. 

•  Case  B  denotes  a  reduction  in  the  parameter  uncertainty  to 
5%  while  the  variability  on  the  power  density  is  sought  to 
be  kept  within  0.02  W  cm-2  as  in  Case  A.  Corresponding  to 
these  specifications,  the  following  parameter  combinations 
and  maximum  realizable  power  densities  are  obtained 
from  Figs.  1 1  and  12:  T  =  100  °C,  pd-  10  atm,  pc  =  15  atm, 
M/>d,max  =  0.65  W  cm-2  (Fig.  1 1(c));  T  =  100  °C,  pd  =  6  atm, 
pc  =  15  atm,  jipd,max  =  0.68  W  cm-2  (Fig.  12(b))  and 
T  =  100°C,/?a  =  10  atm, pc  =  15  atm,  /xpd,max  =  0.65  W  cm-2 
(Fig.  12(c)).  The  operating  parameters  that  lead  to  the  max¬ 
imum  mean  power  density  of  0.68  W  cm-2,  noted  above 
from  Fig.  12(b),  therefore,  constitute  the  robust  parameter 
design  for  Case  B.  For  this  design,  and  the  variability  of 
0.02  W  cm-2,  the  coefficient  of  variance  on  the  power  density 
is  2.9%. 

•  Case  C  explores  the  impact  of  accommodating  a  greater  vari¬ 
ability  in  the  power  density,  while  the  parameter  uncertainty 
is  maintained  at  5%  as  in  Case  B.  The  reader  may  verify  that 
the  robust  operating  conditions  that  lead  to  maximum  mean 
power  density  correspond  to  the  parameter  combination  in 
Fig.  12(b),  as  in  Case  B:  T  =  100°C,pa  =  6atm,pc  =  15  atm,  for 
which  the  maximum  realizable  power  density  is  0.90  W  cm-2, 
and  its  coefficient  of  variance  is  5.6%. 

A  comparison  of  the  three  design  cases  indicates  that  reducing 
the  parameter  uncertainty  from  10%  (Case  A)  to  5%  (Case  B) 
leads  to  a  70%  increase  in  the  maximum  realizable  power  density 
from  0.40  to  0.68  W cm-2.  Further,  by  allowing  for  a  greater 
variability  of  0.05  W  cm-2  in  the  power  density  (Case  C)  while 
keeping  the  parameter  uncertainty  at  5  % ,  the  mean  power  density 
is  increased  to  0.90  W  cm-2 — an  increase  of  32%  relative  to  the 
0.68  W  cm-2  in  Case  B.  Of  the  three  cases,  the  design  for  Case  B 
represents  the  least  coefficient  of  variance  of  the  power  density. 
The  contour  maps  in  Figs.  11  and  12  thus  provide  for  ready 
exploration  of  tradeoffs  in  arriving  at  robust  designs  based  on 
application  requirements. 

The  results  presented  here  were  based  on  parametric  studies 
over  a  selected  range  of  three  principal  operating  parameters 
in  order  to  illustrate  the  methodology  of  fuel  cell  design  under 
uncertainty  using  stochastic  simulations.  The  stochastic  analysis 
methodology  may  be  extended  to  design  considerations  involv¬ 
ing  parameters  other  than  the  power  density,  and  moreover, 
other  considerations,  such  as  those  involved  in  the  applica¬ 
tion  of  high  cell  pressures,  may  be  incorporated  as  additional 
objectives  and  constraints  in  the  design  endeavor.  The  stochas¬ 
tic  analysis  methodology  was  demonstrated  in  this  study  using  a 
one-dimensional  single-phase  model  of  a  PEM  fuel  cell,  which 
served  to  illustrate  the  principal  trends  and  results  that  are  appli- 
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Fig.  12.  Contour  maps  representing  the  maximum  mean  value  of  power  density,  in  a  two-dimensional  field  of  the  standard  deviation  of  power  density  and  input 
uncertainty,  for:  (a)-(c)  various  values  of  anode  pressure  and  (d)-(f)  various  values  of  cathode  pressure.  Other  operating  parameters  are  given  in  Table  1 . 


cable  to  practical  fuel  cells.  The  methodology  can  be  readily 
extended  to  more  sophisticated  multi-dimensional  and  multi¬ 
phase  models  of  fuel  cells  and  their  stacks  with  an  associated 
increase  in  the  computational  effort.  These  enhancements  may 
be  considered  in  a  future  study.  In  a  future  work,  the  stochastic 
modeling  framework  may  also  be  interfaced  with  a  numerical 
optimization  scheme  to  provide  a  robust  design  tool  for  stochas¬ 
tic  optimization  under  uncertainty. 

5.  Conclusions 

The  paper  presented  detailed  discussion  on  the  framework 
of  a  sampling-based  stochastic  analysis  for  proton  exchange 
membrane  fuel  cells,  where  some  of  the  operating  and  design 
parameters  (input  parameters)  are  uncertain.  The  results  were 
presented  in  terms  of  mean  and  standard  deviation  of  power 
density,  which  represent  a  performance  measure  and  its  vari¬ 
ability.  The  methodology  and  the  results  provide  for  a  valuable 
tool  for  fuel  cell  design  under  uncertainty.  Example  cases  illus¬ 
trating  the  use  of  the  results  for  robust  design  under  uncertainty 
were  presented.  The  stochastic  analysis  framework  presented  in 
the  paper  may  be  used  for  cases  where  other  input  parameters 
are  considered  to  be  uncertain.  The  stochastic  model  can  also  be 


used  as  a  basis  for  numerical  optimization  under  uncertainty  for 
fuel  cells. 
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